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ABSTRACT 

Aims. We study the long-term time-scale (i.e. period comparable to the orbital period of the outer perturber object) transit timing 
variations in transiting exoplanetary systems which contain a further, more distant (a 2 >> fli) either planetary, or stellar companion. 
Methods. We give an analytical form of the O - C diagram (which describes such TTV-s) in trigonometric series, valid for arbitrary 
mutual inclinations, up to the sixth order in the inner eccentricity. 

Results. We show that the dependence of the O - C on the orbital and physical parameters can be separated into three parts. Two 
of these are independent of the real physical parameters (i.e. masses, separations, periods) of a concrete system, and depend only on 
dimensionless orbital elements, and so, can be analysed in general. We find, that for a specific transiting system, where eccentricity 
(ei) and the observable argument of periastron (a>i) are known e.g. from spectroscopy, the main characteristics of any, caused by a 
possible third-body, transit timing variations can be mapped simply. Moreover, as the physical attributes of a given system occur only 
as scaling parameters, the real amplitude of the O - C can also be estimated for a given system, simply as a function of the m 3 /P 2 
ratio. We analyse the above-mentioned dimensionless amplitudes for different arbitrary initial parameters, as well as for two particular 
systems CoRoT-9b and HD 80606b. We find in general, that while the shape of the O - C strongly varies with the angular orbital 
elements, the net amplitude (departing from some specific configurations) depends only weakly on these elements, but strongly on the 
eccentricities. As an application, we illustrate how the formulae work for the weakly eccentric CoRoT-9b, and the highly eccentric HD 
80606b. We consider also the question of detection, as well as the correct identification of such perturbations. Finally, we illustrate 
the operation and effectiveness of Kozai cycles with tidal friction (KCTF) in the case of HD 80606b. 

Key words, methods: analytical - methods: numerical - planetary systems - binaries: close - Planets and satellites: individual: 
CoRoT-9b - Planets and satellites: individual: HD 80606b 



' 1. Introduction mid-transit times with respect to the cycle numbers we get the 

O - C diagram which has been the main tool for period studies 
The rapidly increasing number of exoplanetary systems, as well by variable sta r observers (not only for eclipsing binaries) for 
. as the lengthening time interval of the observations naturally more than a cen tury. Consequently, the effect of the various types 
leads to the search for perturbations in the motion of the known of period variations (both rea i and apparent) for the O - C dia- 
planets which can provide the possibility to detect further plan- gram were a]ready widely studied - m the last one hundred years, 
etary (or stellar) components in a given system, and/or can pro- Some of these are less relevant m the case of transiting exo- 
duce further information about the oblateness of the host star (or planets> but others ^ important . For examp i e , the two classical 
the planet), or might even refer to evolutionary effects. cases are the simple geometrical light-time effect (LITE) (due 
The detection and the interpretation of such perturbations in t0 a f urt her, distant companion), and the apsidal motion effect 
the orbital revolution of the exoplanets usually depends on such (AME) (due to both the stellar oblateness in eccentric binaries 
methods and theoretical formulae which are well-known and and tne relativistic effect). Due to its small amplitude however, 
have been applied for a long time in the field of the close eclips- LITE, which has been widely used for identification further stel- 
ing and spectroscopic binaries. We mainly refer to the methods i ar companions of many variable stars (not only in case of eclips- 
developed in connection with the observed period variations in i ng binaries) since the papers of Chandler (1892); Hertzsprung 
eclipsing binaries. These period variations manifests themselve (1922); Woltjer (1922); Irwin (1952), is less significant in the 
in the departure of the occurrence of an eclipse event (either tran- C ase of a planetary-mass wide component. Nevertheless in the 
sit, or occultation) from its predicted time. Applying the nomen- rece nt years there have been some efforts to discover exoplan- 
clature of exoplanet studies this phenomenon is called transit e ts on this manner (see e.g. Silvotti et al. 2007, V391 Pegasi, 
timing variation (TTV). Plotting the observed minus calculated Deeg et al. 2008, CM Draconis, Qian et al. 2009, QS Virginis, 
Lee et al. 2009, HW Virginis). The importance of AME in close 
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exoplanet systems has been investigated in several papers (see 
e.g. Miralda-Escude 2002; Heyl & Gladman 2007; Jordan & 
Bakos 2008, and further references therein). Note, that this latter 
effect has been studied since the pioneering works of Cowling 
(1938); Sterne (1939), and besides its evident importance in the 
checking of the general relativity theory, it plays also an impor- 
tant role in studying the inner mass distributions of stars (via 
their quadrupol moment), i.e. it provides observational verifica- 
tion of stellar modells. 

Nevertheless, besides their similarities, there are also sev- 
eral differences between the period variations of close binary 
and multiple stars, and planetary systems. For example, most of 
the known stellar multiple systems form hierarchic subsystems 
(see e.g. Tokovinin 1997). This is mainly the consequence of the 
dynamical stability of such configurations, or put another way, 
the dynamical instability of the non-hierarchical stellar multiple 
systems. In contrast, a multiple planetary system may form a 
stable (or at least long-time quasi-stable) non-hierarchical con- 
figuration as we can see in our solar system, and furthermore, 
the same result was shown by numeric integrations for several 
known exoplanet systems (see e. g Sandor et al. 2007). Due to 
this fact we can expect several different configurations, the ex- 
amination of wich has, until now not been considered in the field 
of the period variations of multiple stellar systems. Examples of 
these are as follows, the perturbations of an inner planet, as well 
as of a companion on a resonant orbit (Agol et al. 2005), or as 
a special case of the latter, the possibility of Trojan exoplanets 
(Schwarz et al. 2009). Recently, the detectability of exomoons 
has also been studied (Simon et al. 2007; Kipping 2009a,b). 

Furthermore, due to the enhanced activities related to ex- 
trasolar planetary searches, which led to missions like CoRoT 
and Kepler which produce long-term, extraordinarily accurate 
data, we can expect in the close future such kinds of observa- 
tions which give the possibility of detecting and studying further 
phenomena which was never observed earlier. For example, it 
is well-known, that neither the close binary stars, nor the hot- 
Jupiter-type exoplanets could have been formed in their present 
positions. Different orbital shrinking mechanisms are described 
in the literature. Instead of listing them, we refer to the short 
summaries by Tokovinin et al. (2006); Tokovinin (2008) and 
Fabrycky & Tremaine (2007). Here we note only, that one of the 
most preferred theories for the formation of close binary stars, 
which also might have produced at least a portion of hot-Jupiters, 
as well, is the combination of the Kozai cycles with tidal friction 
(KCTF). The Kozai resonance (recently frequently referred to 
as Kozai cycle[s]) was first described by Kozai (1962) investi- 
gating secular perturbations of asteroids. The first (theoretical) 
investigation of this phenomenon with respect to multiple stellar 
systems can be found in the studies by Harrington (1968, 1969); 
Mazeh & Saham (1979); Soderhjelm (1982). A higher, third or- 
der theory of Kozai cycles was given by Ford et al. (2000), while 
the first application of KCTF to explain the present configura- 
tion of a close, hierarchical triple system (the emblematic Algol, 
itself) was presented by Kiseleva et al. (1998). According to this 
theory the close binaries (as well as hot-Jupiter systems) should 
have originally formed as significantly wider primordial bina- 
ries having a distant, inclined third companion. Due to the third 
object induced Kozai oscillation, the inner eccentricity becomes 
cyclically so large that around the periastron passages, the two 
stars (or the host star and its planet) approach each other so 
closely that tidal friction may be effective, which, during one 
or more Kozai cycles, shrinks the orbit remarkably. Due to the 
smaller separation, the tidal forces remain effective on a larger 
and larger portion of the whole revolution, and finally, they will 



switch off the Kozai cycles, producing a highly eccentric, moder- 
ately inclined, small-separation intermediate orbit. After the last 
Kozai cycle, some additional tidally forced circularization may 
then form close systems in their recent configurations. 

Nevertheless, independently from the question, as to which 
mechanism(s) is (are) the really effective ones, up to now we 
were not able to study these mechanisms in operation, only the 
end-results were observed. This is mainly the consequence of se- 
lection effects. In order to study these phenomena when they are 
effective, one should observe the variation of the orbital elements 
of such extrasolar planets, as well as binary systems which are 
in the period range of months to years. The easisest way to carry 
out such observations is the monitoring of the transit timing vari- 
ations of these systems. However, there are only a few known 
transiting extrasolar planets in this period regime. Furthermore, 
although binary stars are also known with such separations, sim- 
ilarly, they are not appropriate subjects for this investigation be- 
cause of their non-eclipsing nature. 

The continuous long-term monitoring of several hundreds of 
stars with the CoRoT and Kepler satellites, as well as the long- 
term systematic terrestrial surveys provide an excellent opportu- 
nity to discover transiting exoplanets (or as by-products: eclips- 
ing binaries) with the period of months. Then continuous, long- 
term transit monitoring of such systems (combining the data 
with spectroscopy) may allow the tracing of dynamical evolu- 
tionary effects (i. e. orbital shrinking) already on the timescale 
of a few decades. Furthermore, the larger the characteristic size 
of a multiple planetary, stellar (or mixed) system the greater the 
amplitude of even the shorter period perturbations in the tran- 
sit timing variations, as was shown in detail in the discussion of 
Borkovits et al. (2003). 

In the last few years several papers have been published on 
transit timing variations, both from theoretical aspects (e.g. (non- 
complete) Agol et al. 2005; Holman & Murray 2005; Nesvorny 
& Beauge 2010; Holman 2010; Cabrera 2010; Fabrycky 2010, 
and see further references therein), and large numbers of papers 
on observational aspects for individual transiting exoplanetary 
systems. Nevertheless, most (but not all) of the theoretical pa- 
pers above, mainly concentrate simply on the detectibility of fur- 
ther companions (especially super-Earths) from the transit tim- 
ing variations. 

In this paper we consider this question in greater detail. We 
calculate the analytical form of the long-period 1 (i. e. with a pe- 
riod on the order of the orbital period of the ternary component 
P2) time-scale perturbations of the O - C diagram for hierarchi- 
cal (i.e. P2 » P\) triple systems. (Note, as we mainly concen- 

1 In this paper we follow the original classification of Brown 
(1936) who divided the periodic perturbations occuring in hierar- 
chical systems into the following three groups: 

- Short period perturbations. The typical period is equal to the 
orbital period Pi of the close pair, while the amplitude is of the 
order (Pi/P 2 ) 2 ; 

- Long period perturbations. This group has a typical period of 
P 2 , and magnitude of the order (P1/P2); 

- Apse-node terms. In this group the typical period is anout 
P\/Pi, and the order of the amplitude reaches unity. 

This classification differs from what is used in the classical plane- 
tary perturbation theories. There, the first two groups termed to- 
gether "short period" perturbations, while the "apse-node terms" 
are called "long period ones". Nevertheless, in the hierarchical sce- 
nario, the first two groups differ from each other both in period and 
amplitude, consequently, we feel this nomenclature more appropri- 
ate in the present situation. 
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trated on transiting systems with a period of weeks to months, 
we omitted the possible tidal forces. However, our formulae can 
be practically applied even for the closest exoplanetary systems, 
because the tidal perturbations become effective usually on a no- 
tably longer time-scale.) This work is a continuation and exten- 
sion of the previous paper of Borkovits et al. (2003). In that paper 
we formulated the long-period perturbations of an (arbitrarily ec- 
centric and inclined) distant companion to the O - C diagram for 
a circular inner orbit. (Note, that our formula is a generalized 
variation (in the relative inclination) of the one of Agol et al. 
2005. For the coplanar case the two results become identical.) 
Now we extend the results to the case of an eccentric inner bi- 
nary (formed either a host-star with its planet, or two stars). As 
it will be shown, our formulae have a satisfactory accuracy even 
for a high eccentricity, such as e\ = 0.9. Note, that in the pe- 
riod regime of a few months the tidal forces are ineffective, so 
we expect eccentric orbits. This is especially valid in the case of 
the predecessor systems of hot-Jupiters, in which case the above 
mentioned theory predicts very high eccentricities. 

In the next section we give a very brief summary of our 
calculations. (A somewhat detailed description can be found in 
Borkovits et al. 2003, 2007, nevertheless for self-consistency of 
the present paper we provide here a brief overview .) In Sect. 
3. we discuss our results, while in Sect. 4 we illustrate the re- 
sults with both analytical and numerical calculations on two in- 
dividual systems, CoRoT-9b and HD 80606b. Finally in Sect. 
5 we conclude our results, and, furthermore, we compare 
our method and results with that of Nesvorny & Morbidelli 
(2008); Nesvorny (2009). 

2. Analytical investigations 

2.1. General considerations and equations of the problem 

As is well-known, at the moment of the mid-transit (which in 
case of an eccentric orbit usually does not coincide with the half- 
time of the whole transit event) 



+ - + 2kn, 
2 



(1) 



where u is the true longitude measured from the intersection 
of the orbital plane and the plane of sky, and k is an integer. 
(Note, since, traditionally the positive z-axis is directed away 
from the observer, the primary transit occurs at u — -n/2.) An 
exact equality is valid only if the binary has a circular orbit, or 
if the orbit is seen edge-on exactly. (The correct inclination de- 
pendence of the occurrence of the mid-eclipses can be found in 
Gimenez & Garcia-Pelayo 1983.) Nevertheless, the observable 
inclination of a potentially month-long period transiting extra- 
solar planet should be close to ; = 90°, if the latter condition is 
to be satisfied. Due to its key role in the occurrence of the tran- 
sits, instead of the usual variables, we use u as our independent 
time-like variable. It is known from the textbooks of celestial 
mechanics, that 



u — —z - Q cos ;', 
P\ 

,l/2„-3/2. 



// 1/Z a- J/Z (l -e 2 r 3/2 (l +ecosv) z -Qcosi, 



(2) 



consequently, the moment of the A^-th primary minimum (tran- 
sit) after an epoch to can be calculated as 



Jt J-n/2 



.2Nn-7i/2 a 3/2 



(l_ e 2)3/2 



du 



H 1 ' 2 [1 + e cos( M -w)] 2 ! _dQ C0S /' 



,3/2 



(l_ e 2 )3 /2 



Pi 



l + -!-ncosj \du. (3) 
fi 1 ' 2 [1 + ecos(t< - co)] 2 \ c\ 



In the equations above c\ denotes the specific angular momen- 
tum of the inner binary, p t is the radius vector of the planet 
with respect to its host-star, while the orbital elements have their 
usual meanings. Furthermore, in Eq. (3) we applied that the true 
anomaly can be written as v = u - co. In order to evaluate Eq. (3) 
first we have to express the perturbations in the orbital elements 
with respect to u. Assuming that the orbital elements (except of 
u) are constant, the first term of the right hand side yields the 
following closed solution 



PiJ 



P 
2^ 



2arctan 



1 - e + cos co 
1 + e 1 + sin co 



± (1 - e z ) 



2n1/2_ 



ecosw 



1 + e sin co 



(4) 



for the two types of minima, respectively. (Here P denotes the 
anomalistic or Keplerian period which is considered to be con- 
stant.) Note, that instead of the exact forms above, widely used 
is the expansion (as in this paper), which, up to the fifth order in 
e, is as follows: 



P m = P S E+ — 
2n 



1 „ I 3 2 1 4 1 ■ 
+ -7r ± 2e cos co + —e + —e I sin2w 

2 \ 4 8 



1 3 1 5 \ , 5 4 . 3 5 

-e + —e cos3w e sin4w± — e cos5<jj 

3 8/ 



32 



40 



(5) 



where P s is the sidereal (or eclipsing) period of, for example, 
the first cycle, and E is the cycle-number. In the present case 
the orbital elements cease to remain constant. Nevertheless, 
as it can be seen from Eqs. (8)-(14), their variation on P 2 
timescale is related to Pi/Pi « 1, which allows the lin- 
earization of the problem, i.e. in such a case Eq. (5) is for- 
mally valid in the same form, but e, co and P s are no longer 
constant. Then a further integration of Eq. (5) with respect 
to v 2 gives the analytical form of the perturbed O - C on P 2 
time-scale. To this, as a next step, we have to calculate the long- 
period and apse-node perturbations in u. Some of these arise 
simply from the similar perturbations of the orbital elements (or 
directly of the e cos co, e sin co functions), while others (we will 
refer to them as direct perturbations in u) come from the varia- 
tions of the mean motion (for more details see Borkovits et al. 
2007) 2 . (In other words this means that in such a case P will no 
longer be constant.) 

At this point, to avoid any confusion, we emphasize that dur- 
ing our calculations we use different sets of the angular orbital 
elements. As we are interested in such a phenomenon which pri- 
marily depends on the relative positions of the orbiting celestial 
bodies with respect to the observer, the angular elements (i.e. u, 
co, i, Q.) should be expressed in the "observational" frame of ref- 
erence having the plane of the sky as the fundamental plane, and 
u, as well as co is measured from the intersection of the binary's 
orbital plane with that plane, while Q is measured along the 
plane of the sky from an arbitrary origin. On the other hand, the 
physical variations of the motion of the bodies depend on their 
relative positions to each other, and consequently, it is more ben- 
eficial (and convenient) to express the equations of perturbations 
of the orbital elements in a different frame of reference, (we shall 
refer to it as the dynamical frame) which depends on the relative 



2 The minus signs on the rhs of Eqs. (10) and (1 1) of Borkovits et al. 
(2007) should be replaced by plus. 
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Fig. 1. The spatial configuration of the system. 



positions of the bodies, independently from any outer observer. 
The fundamental plane of this latter frame of reference is the in- 
variable plane of the triple system, i.e. the plane perpendicular to 
the net angular momentum vector of the complete triple system. 
In this frame of reference the longitude of the ascending node ( h) 
gives the arc between the sky and the corresponding intersection 
of the two orbital planes measured along the invariable plane, 
while the true longitude (w) and the argument of periastron (g) 
is measured from that ascending node, along the respected or- 
bital plane. In order to avoid a further confusion, the relative 
inclination of the orbits to the invariable plane is denoted by j. 
The meaning and the relation between the different elements can 
be seen in Fig. 1, and listed also in Appendix A. 



2.2. Long period perturbations 

From now on we refer to the orbital elements of the inner binary 
(i.e. the pair formed by either a host star and the inner planet, 
or two stellar mass objects) by subscript i, while those of the 
wider binary (i.e. the orbit of the third component around the 
centre of mass of the inner system) by subscript 2- The differ- 
ential perturbation equations of the orbital elements are listed in 
Borkovits et al. (2007) 3 . (In that paper, from practical consid- 
erations, we did not restrict ourselves to the most usual rep- 
resentation of the perturbation equations, i.e. the Lagrange 
equations with the perturbing function, rather we used the 
somewhat more general form, expressing the perturbations 
with the three orthogonal components of the perturbing 
force. Nevertheless, as far as only conservative, three-body 
terms are considered, the two representations are perfectly 
equivalent.) In order to get the long-period terms of the pertur- 
bation equations, the usual method involves averaging the equa- 
tions for the short-period (» Pi) variables, which is usually the 
mean anomaly (/i) or the true anomaly (vi) of the inner binary, 
but in our special case it is the true longitude u\ . This means that 
we get the variation of the orbital elements by averaged for an 
eclipsing period. Furthermore, in the case of the averaged equa- 
tions we change the independent variable from u i to (the aver- 



3 Note, in the denominator of the Eqs. (14) and (15) of that paper 
for e and lj a closing bracket is evidently missing, and furthermore, the 
equation for Q should be divided by e for the correct result. 



aged) V2, by the use of 



du\ 
dv 2 



C2 p\ 



1 Qi cos i\ 



(6) 



from which, after averaging we get: 



dvl 



p\ 



dQ, 

dV*, 



3/2 I ~ 

"i ^202(1 -e\) 

p 2 a- e jy/2 dQi 

Pi (1 +e2COSV2) 2 dv2 



■ cos l\ 



COS Z]. 



(7) 



In the following we omit the overlining. Furthermore, as one 
can see from Eq. (13) below, the second term in the r.h.s of (7) 
is of the order P1/P2, and consequently can be neglected. So for 
the long-period perturbations of the orbital elements of the close 
orbit we get as follows: 



doi 

dv'2 



^ = 0, 



dv2 



dv 2 



dhi 
dv 2 



dji 
dv2 



(8) 



= A L (l-^) I/2 ei (l+ e2 C0SV2) 



[(l-/ 2 )sin2gi 



--(1 +/) 2 sin(2v 2 + 2# 2 - 2g0 
+ i(l -/) 2 sin(2v 2 + 2g 2 + 2g I ) 



= A L (l-e 2 ) 1/2 (l+e 2 cosv 2 ) 



(9) 



■(I-/ 2 ) 



cos 2gi + - cos(2v 2 + 2g 2 ) 



+ -(l+I) 2 cos(2v 2 + 2g2-2 gl ) 
+ i(l -/) 2 cos(2v 2 + 2g 2 + 2gi) 



dhj_ 

dv 2 



cosji, 



(10) 



-AUl-eh- 1 ' 2 



sin i„ 



sin;j 



-(l + <?2 COS V2) 



X 



-\l + -e 



-ef/cos 2gi 

-^l + ^ 2 J/cos(2v 2 + 2g 2 ) 
+ X -e\{\ +/)cos(2v 2 +2g2-2gi) 



--e\{\ -/)cos(2v 2 + 2g 2 + 2gi ) 



= A L (l-e 2 ) l/2 sin/ m (l +e 2 cosv 2 ) 



-|l + -^|sin(2v 2 + 2g 2 ) 
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dV2 



dQ.i 

dv 2 



dV2 



+e 2 7sin2gi 

--e\{\ +/)sin(2v2 + 2 g2 -2 gl ) 

1 9 

+ -e 2 (l - /) sin(2v 2 + 2g 2 + 2 gi ) 
A L (1 -e\) ll2 {\ + e 2 cosv 2 ) 



x -i - I / z - ^ 



v 5 \ 3 
+ (l-7 2 ) 



cos 2gi + - cos(2v 2 + 2g 2 ) 



1 



+-(l+/) z cos(2v 2 + 2g 2 -2gi) 
+ -/) 2 cos(2v 2 + 2^2 + 2^i) 



dOj 

dv 2 



cos Ii, 



=A L (1-^) 



2n-1/2^ 



sin;i 



-(1 +e2COSV2j 



-II + -e z |/cos(wi -gi) 

+ + 5 e ^ (: /,COS(2 '' 2 + '' a:: ' ' 1 
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-e^cosCwi + gi) 
1 2 

+ -ef(l + 7)cos(2v 2 + 2g 2 -(Oi- gi) 
1 2 

--ef(l -7)cos(2v 2 + 2^2 + ^1 + gi) 



A L (1 - e 2 ) 1/2 sin / m (l + e 2 cos v 2 ) 
2(. 3 

+^|1 + i - /)sin(2r 2 + 2^ - w, +,1;,) 



~5 \ + 2 ei l /sin( ^ Wl 



+ 7 I 1 + 2 e i) (1 + 7 ) sin ( 2v 2 + 2g 2 + wi - gi) 



+e 7sin(<wi + gi) 

1 , 

--e^l + 7) sin(2v 2 + 2^2 - wi -gi) 
1 9 

+ -ef (1 - 7) sin(2v 2 + 2g 2 + wi + gi) 



and, finally, the direct term is as follows: 



dV2 



A L (l-e 2 ) 1/2 (l+ e2 cosv2) 



+ ^( 1 " /2 ) e i /2(ei)cos2 ^ 1 



(11) 



(13) 



(14) 



+ -(l-7 2 )/i( ei )cos(2v 2 +2g2) 
+ ^(1 + 7) 2 e 2 / 2 ( ei ) cos(2v 2 + 2g 2 - 2 gl ) 



40 

+ ^(1 - 7) 2 e 2 / 2 ( ei ) cos(2v 2 + 2g 2 + 2^)} , 



where 



15 P± 

8m 123 /V 2 ' 

and 

f ,, , 25 2 15 4 95 6 

. 31 2 23 4 
/2(e)= 1 + 5i e + 48 e ' 



(15) 
(16) 

(17) 
(18) 



Furthermore, i m denotes the mutual inclination of the two orbital 
(12) planes, while 



7 = cos z' m , 



(19) 



and m 12 3 stands for the total mass of the system. Note, formal in- 
tegration of the first five of the equations above (i.e. those which 
refer to the orbital elements in the dynamical system, Eqs. [8]- 
[11]) reproduce the results of Soderhjelm (1982). 

Strictly speaking, some of the equations above are valid only 
in that case when both orbits are non-circular, and the orbital 
planes are inclined to each other, and/or to the plane of the sky. 
For example, if the outer orbit is circular, neither v 2 , nor g 2 has 
any meaning. Nevertheless, their sum i. e. v 2 + g 2 is meaning- 
ful as before. Furthermore, although in the case of circular in- 
ner orbits the derivative of g\ has no meaning, yet the derivative 
of e\ cosgi, e\ singi (or e\ cos wi, e\ sinwi), i. e. the so-called 
Lagrangian elements can be calculated correctly. (Note, instead 
of the "pure" derivatives of e\ and g\ (or a>\) these latter occur 
directly in the O - C.) Similar redefenitions can be done in the 
case of coplanarity. So, for practical reasons, and for the sake of 
clarity, we retain the original formulations even in such cases, 
when it is formally not valid. 

As one can see, there are some terms on the r.h.s. of these 
equations which do not depend on v 2 . Primarily, these terms 
give the so-called apse-node time-scale contribution to the vari- 
ation of the orbital elements and the transit timing variations. 
Nevertheless, in order to get a correct result for the long-term 
behaviour of the orbital elements these terms must neverthe- 
less be retained in our long-term formulae. Note, these terms 
were calculated for tidally distorted triplets in Borkovits et al. 
(2007). Such formulae (after an omission of the tidal terms) are 
also valid for the present case in the low mutual inclination (i.e. 
approx. 7 2 > |) domain. (Formulae valid for arbitrary mutual 
inclinations will be presented in a subsequent paper.) 

Carrying out the integrations, all the orbital elements on the 
r.h.s of these equations with the exception of v 2 are considered 
as constants. This can be justified for two reasons. First, as one 
can see, for P 2 » Pi the amplitudes of the long-period per- 
turbative terms (A L ) remain small (which is especially valid for 
the case where the host star is orbited by two planets, when 
mi « OT123) 4 , and second, although the amplitude of the apse- 
node perturbative terms can reach unity, the period is usually so 



20 



4 This assumption is analogous to that of the classical low or- 
der perturbation theory where the squares of the orbital ele- 
ment changes are neglected, as they should be proportional to 
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long (~ P\jP\, see e.g. Brown 1936), that its contribution can 
be safely ignored during one revolution of the outer object. 5 The 
final result of such an analysis is an analytical form of the tran- 
sit timing variations, i.e. the long-period O - C diagram. Here 
we give the result up to the first order in the inner eccentricity, 
while a more extended result up to the sixth order in the inner 
eccentricity can be found in appendix B, where we give also the 
perturbation equations directly for the e'" cos nwi,e™ sin nuj\ ex- 
pressions. 



0-Cp 



Pi 
Yn 



4 5 2 6 . 
- + -e, + -e\ sin wi 



) 



[l 2 - 1 - ] jM+ 1 -(l-I 2 )S(2v 2 + 2g2) 

^e 2 cos 2g t + 2ei sin(wi - 2g x ) 

{\-I 2 )M+ l -(\+I 2 )S{2v 2 + 2g 2 ) 
51 , 

— e\ sin2gi + 2ei cosfwj - 2 gl ) IC(2v 2 + 2g 2 ) 



+ cot;'i sin; m {- - (1 + 2e\ sinwi)cosM m i/ 



M - -S(2v 2 + 2g 2 ) 



+ -(1 + 2e\ sinw 1 )sinM ml C(2v2 + 2g 2 ) 



mi a 2 smi 2 



(l - e^sin^ + u> 2 ) 



m m c 



1 + e 2 cos v 2 



+0(ei), 



(20) 



where 
M 



1 + e 2 cos v 2 dv 2 

— v 2 — l 2 + e 2 sin v 2 

3 1 

= 3e 2 sin v 2 - -e\ sin 2v 2 + -e\ sin 3v 2 + 0(e 2 ) 



(21) 



= 3e 2 1 



3 \ 9 53 

- -e 2 sin l 2 + —e\ sin 2l 2 + —el sin 3l 2 + <9(e 2 ), 
o J 4 24 



(22) 



and 



<S(2v 2 + x) = sin(2v 2 + x) + e 2 sin(v 2 + x) + -e 2 sin(3v 2 + x) 



(l - 4e 2 ) sin(2/ 2 + x) 



('Mperturber/'Mstar) 2 , but with the difference that here the small param- 
eter is the ratio of the semi-major axes aja 2 instead of the mass 
ratio. Consequently, this assumption remains valid even if the per- 
turber would be a sufficiently distant stellar-mass object. 

3 Note, this is not necessarily true in the presence of other perturba- 
tive effects. Nevertheless, in such cases one can assume, that if there is 
a physical effect producing apse-node time-scale perturbations having 
a period comparable to the period of the long-period perturbations aris- 
ing from some other sources, then the amplitudes of these long-period 
perturbations are usually so small with respect to the amplitudes of the 
apse-node perturbations that their effect can be neglected. 



+e 2 
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C(2v 2 + x) = cos(2v 2 + x) + e 2 cos(v 2 + x) + -e 2 cos(3v 2 + x) 

= (l -4£>2)cos(2Z 2 +x) 
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4n 
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, I () e l I cos(3/ 2 + v) 
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+e 



1 17 
-- cos x + — cos(4/ 2 + x) 
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cos(/ 2 - x) H cos(5Z 2 + x) 
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furthermore, 



"ml = Wl -gi, 



(24) 



+ 0{e\), 



(25) 



i.e. u m \ is the angular distance of the intersection of the two or- 
bits from the plane of the sky, or, in other words, the longitude 
of the (dynamical) ascending node of the inner orbit along the 
orbital plane, measured from the sky. 

Note, the superscripts refer to the exoplanetary transits (pri- 
mary minima), and the subscripts to the secondary occultations 
(secondary minima). Furthermore, we assumed the formally 
second-order ^e\, and ^e\ terms, to be first order, as their val- 
ues exceed e\ for medium eccentricities. We included also the 
pure geometrical light-time contribution in the last row. Here c 
denotes the speed of light. The minus sign arises because this 
term reflects the motion of the inner pair around the common 
centre of mass, whose true longitude differs by n from the one of 
the third component. In Eq. (21) the mean anomaly of the outer 
body (l 2 ) appears because of (the constant part of) the difference 
between the anomalistic P and sideral (eclipsing, or transiting) 
P s period included into the first term in the r.h.s. of Eq. (5). 



3. Discussion of the results 

As one can easily see, the result for a circular inner orbit (i.e. 
e\ — 0) is identical with the formula (46) of Borkovits et al. 
(2003). Consequently, the discussion given in Sect. 4 of that 
paper is also valid. Nevertheless, as one can see in Figs. 2, 3, 
a significant inner eccentricity produces notably higher ampli- 
tudes, and consequently, is easier to detect. Furthermore, some 
other attributes of the transit timing variations also change dras- 
tically. For example, in contrast to the previously studied copla- 
nar (/ = ±1), circular inner orbit (e\ =0) case, (Borkovits et al. 
2003; Agol et al. 2005), as long as the inner orbit is eccentric, 
the dynamical term does not disappear even if the outer orbit is 
circular (e 2 = 0). A further important feature for eccentric inner 
orbit, in that the amplitude, the phase and the shape of the O - C 
variations cease to be depend simply on the physical (i.e. rela- 
tive) positions of the celestial bodies, but also on the orientation 
of the orbit with respect to the observer. These latter elements 
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(i. e. e\ sinwi, e\ cos u>\ and their combinations) can also be de- 
termined from radial velocity measurements, as well as from the 
shape of the transit light-curves, and, in the case of the possible 
detection of the occultation, or (in the case of eclipsing binaries), 
of the secondary minimum, from the time delay between the two 
different eclipsing events. So, while the relative, i.e. physical an- 
gular parameters (i.e. periastron distances from the intersection 
of the two orbits, g\, g2, and mutual inclination z' m ) cannot be 
aquired from other, generally used methods (e.g. light-curve, or 
simple radial velocity curve analysis), these observational geo- 
metrical parameters could be determined from other sources of 
information, and then simply can be built into such a fitting algo- 
rithm, which was described in Borkovits et al. (2003), or could 
be included into such procedures which were presented by Pal 
(2010). 

In the following, as we are mainly interested in transiting 
systems, where i\ as 90° (which is especially true for relatively 
longer period, i.e. distant transiting exoplanets), we omit terms 
multiplied by cotz'i, i.e. terms arising from nodal motion. 

In Borkovits et al. (2003) we considered only triple stellar 
systems, where the three masses were usually nearly equal, and 
the inner period was of the order of a few days. In such cases 
light- time term dominates. The amplitude of the light-time effect 
is simply 



m3 fl2Sin;2 
m m c 
mi [Gm 



(l -e^cos 2 ^) 



1/2 



(Gm m \ l/3 sini 2 2/3/, 2 2 \ 

-4™3 „ ; „, D 2/3 A e 2 cos 2^ 



1/2 



«123 



1.1 x lO- 4 -^sim 2 P2 /3 (l 



m 



2/3 

123 



(26) 



where masses should be expressed in solar units, and P 2 in days. 

The amplitude of the dynamically forced O - C, and con- 
sequently, the detectability limit of such perturbations depends 
on almost the all dynamical as well as geometrical variables. 
So, we can give only some limits on the detectability limit. 
Nevertheless, for an easier, and somewhat general study, we sep- 
arate the physical and geometrical variables from each other, 
and furthermore, we separate also the elements of the inner orbit 
from those of the outer perturber (with the exception of the mu- 
tual inclination). In order to do this, we introduce the following 
quantities, all of which depend on eccentricity (ei), the two types 
of periastron arguments (gi, wi), and mutual inclination (; m ), via 
its cosine: 



'2 5 

" = * 5 + 4 e i ^ +cos2 -?i) + e i 



■ sinw! + sin(w! - 2gi) 
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4' 



cos2gi) 
sinwi - sin(wi - 2g\) 
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5 , 

—e l sin2gi + e\ cos(wi - 2g\) 



(27) 
(28) 



4 5 



1 
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■2e x 



sin(wi - 2g\) - - sin co\ 



[4 5 2 
\5 + 2 e ' 



(1 - cos 2gi) 



-2e, 



— sin cji - sin(o>i - 2g\) 



}} ('-;)'" 



(29) 



Then 

O - Cdyn 

where 



2n 



yM +^a 2 +0 1 S{2v 2 + < 



(f> = 2g2 + arctan 



t 
a 



or, in trigonometric form 

O - C dyn ~ ^A* L ^ A„ sin(nv 2 + 4>n), 



(30) 



(31) 



(32) 



where (up to third order in e 2 ) 

Ai = ( 1 - 4) a a/ A | + A m + 2A s A M cos cf>, (33) 
A 2 = (l -4)~ V2 + ^4 a m + \ e l A sA M cos ~(j>, (34) 
A 3 = (l-ej) 3/2 \e 2 Ja 2 s + ^e 4 2 A 2 M + ^e\A s A M cos 0, (35) 

(36) 



sin <p 

(pi = arctan j— , 

cos0+^ 

sin0 



cp2 = arctan 



03 = arctan 



cos <f>+\e 2 2 Tj' 



sin< 



cos 0+14^' 



and 

A s = yja 2 +j3 2 , 
A M = 3y, 



and, furthermore, 



Al=A L (l-ef 2 . 
Note, for a circular outer orbit (e 2 = 0) 



A3 
A 2 



0, 
As. 



(37) 
(38) 

(39) 
(40) 

(41) 



(42) 
(43) 



As one can see, both sets of amplitudes, i.e. Am,s and A 1,2,3 
are independent of the real physical parameters of the current 
exoplanetary system. Furthermore, the dependence from the el- 
ements of the outer orbit (i.e. e2, g 2 ) appear only in the A 1,2,3 
Fourier amplitudes. The masses (or more accurately mass ra- 
tios), and the periods and period ratios (and indirectly the phys- 
ical size of the system) occur only as a scaling parameter. In 
this way, the following general statements are valid for every 
hierarchical triple systems (as far as the initial model assump- 
tions are valid). In order to get the real, i.e. physical values for 
the O - C amplitudes in a given system, one must multiply the 
general system-independent, dimensionless amplitudes with the 

system specific number j^- jr ■ 



16tt P2 m\+ni2+mi ' 

Considering first the (nearly) coplanar case, i.e. when / 
-A, then the (half-)amplitude of the two sinusoidals becomes: 



8/ 2 si/2/ 25 
eiA M o ~ 5 ( ! _e ij 1 + y 
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25 , 

1 sinwi + — e\. 



(44) 



(45) 
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This illustrates the above-mentioned fundamental differences to 
the previously-studied coplanar, e\ = case (Borkovits et al. 
2003; Agol et al. 2005), as according to our new result for ec- 
centric inner orbits, the dynamical term does not disappear even 
if the outer orbit is circular (e 2 = 0). Furthermore, the amplitudes 
strongly depend on the orientation of the orbital axis with respect 
to the observer. When the apsidal line coincides (more or less) 
with the line of sight (i.e. lo\ = ±90°) there can be very signifi- 
cant differences both in the shape and amplitude of the primary 
(transit) and secondary (occultation) O - C curves. However, 
when the apsidal line lies nearly in the sky, then these differences 
disappear. A further interesting feature of the (±>\ = +90° config- 
uration, is that for eclipse events which occur around apastron 
there is a full square under the square-root sign in the .S-term, 
with the root of e\ = 0.8, which means that in this situation this 
term would disappear. Nevertheless, for such a high eccentricity 
the first order approximation is far from being valid, as is illus- 
trated in Fig. 2. 

For highly inclined (/ « 0) orbits the (half-)amplitudes are 
as follows: 



Fig. 2 shows the inner eccentricity (ei) dependence of the 
amplitudes for some specific values of the other parameters. In 
general, one sees that the e\ increase with amplitude, but this 
growth usually remains within one order of magnitude. So, al- 
though the shapes of the individual O - C curves may differ 
significantly, the net amplitudes vary over a narrow range. The 
graphs also suggest another, perhaps suprising fact, that the am- 
plitude of O - Cp2dyn depends only weakly on the mutual incli- 
nation (; m ). This is especially valid for medium inner eccentric- 
ities, since in this case, at least for the cases shown in Fig. 2, 
all Am,s amplitudes have similar values. Moreover, the numeri- 
cally generated sample O-C curves of Fig. 3, also suggest such 
a conclusion. Nevertheless, there are several exceptions. Some 
particular configurations some of the amplitudes, or even both 
of them (and consequently, the corresponding terms) can dis- 
appear. Such a situation will be discussed in Sect. 4. 1 . We will 
return to this question in detail in Sect. 4, where the dependence 
of the amplitudes on other parameters will be also studied for the 
systems of CoRoT-9b and HD 80606b. 

Returning now to the LITE amplitude, and comparing it with 
the dynamical case, an increment of P2 (keeping Pi constant) re- 
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furthermore, the phase of the <S-term is simply 
= 2g2- 



(46) 



(47) 



(48) 



In this case a further parameter, namely, the periastron distance 
of the inner planet from the intersection of the two orbital planes 
(gi) also plays an important role. 

Finally we also mention a very specific case, namely for the 
maximum eccentricity phase of the Kozai mechanism driven 
e-cycles. During this phase, cos2gi takes one definite value, 
namely cos 2gi = -l. Furthermore, the mutual inclination of the 
two orbits here reaches its minimum. The actual value depends 
on both the maximum mutual inclination ; m , or more strictly, 
on j\, and the minimal inner eccentricity e\. Nevertheless, in 
the case of an initially almost circular inner orbit, the mini- 
mum mutual inclination is almost independent of its maximum 
value, and takes j\(*s i m ) « 39°23 (or its retrograde counterpart, 
f m ) x !40°77), i.e. I 2 = 3/5. For this scenario: 



more distant systems the pure geometrical effect tends to exceed 
the dynamical one, as is the case in all but one (A Tau, see e. g. 
Soderhjelm 1975) of the known hierarchical eclipsing triple stel- 
lar systems. Nevertheless, as we will illustrate in this paper (see 
Fig. 5), we have a good chance of finding the opposite situation 
in some of the recently discovered transiting exoplanetary sys- 
tems. We can make the following crude estimation. Consider a 
system with a solar-like host star, and two approximately Jupiter- 
mass companions) (OT123 ~ m\ = 1M , m 2 — m 3 = 1O~ 3 M ) 
choosing A = 10~ 3 day for the case of a certain detection, then 
the LITE term for the most ideal case gives P2 > I0 6 d a 2 700y. 
Alternatively, setting m 3 = 1O~ 2 M , and allowing A = 10~ 4 for 
detection limit, the result is P2 > I0 3 d « 2.7y. Similarly, for 
two Jupiter-mass planet, in the coplanar case, for small e\, the 
A = 10~ 3 day limit gives 



P 2 <2xl0 3 m 3 e 2 (l-eiy V2 P 2 



(51) 



condition for the detectability of a third companion by its long- 
term dynamical perturbations. For the perpendicular case the 
same limit is 



(52) 
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n ± -e\ sino>i ]e 2 , (49) 



/ 2x1/2/ 4 5 2 34 \ 
aKozai.90 = \ l - e i) I 25 _ 3 e i * 25 e ' slnWl )' 



/ ,\l/2 2 / — 

j8 Ko zai,90 = +(1 -e x ) -V15eicoswi. 



(50) 



In order to get a better overview of the parameter depen- 
dence of the formulae above, we investigate the Am,s< as we ll 
as the A 1,2,3 amplitudes graphically. Due to the complex depen- 
dence of these amplitudes on many parameters, it is difficult to 
give general statements. Therefore we investigate only the de- 
pendence of the amplitudes on the inner eccentricity, while the 
effects of other parameters will be considered in Sect. 4 for spe- 
cific systems, where some of the parameters can be fixed. 



P 2 < Ixl0 3 m 3 (l-^)" 3/2 (l + ^2jP?. 

We note that for m 3 «m\ the above equations are linear for m 3 , 
so it is very easy to give the limiting period P 2 in the function 
of m 3 . Nevertheless, we emphasize again that there are so many 
terms with different periodicity and phase, that these equations 
give only a very crude, first estimation. For Pi =1, 10, 100 days, 
(at zero outer eccentricity) gives P2 <0.5, 50, 5 000 days, respec- 
tively, 

These results refer to the total amplitude of the O-C curve, 
i.e. the variation of the transit times during a complete revolution 
of the distant companion, which can take as long as several years 
or decades. Naturally, the perturbations in the transit times could 
be observed within a much shorter period, from the variation of 
the interval between consecutive minima. This estimation can 
be calculated e.g from Eq. (30), or even directly from Eq. (15). 
According to the meaning of the O - C curve, the transiting or 
eclipsing period (P) between two consecutive (let us assume the 
n-th and n + 1-th) minima is 

P„ = t n+i -t„ = (O-C) OWi) -(O-C) (t n ) 
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Fig. 2. Left panels: The inner eccentricity (e{) dependence of Am, A$ amplitudes of M, S functions of long-period dynamical part of O - C for 
specific values of some parameters. Note, to compare the M and <S terms, the former should be multiplied by e 2 . The thin lines (indexed by ' 1') 
refer to the first order approximation, while the thick ones (index '6') to the sixth one. Right panels: The corresponding Ai^i amplitudes of the 
trigonometric representation (Eq. 32) of the O — C for two different outer eccentricities (ei = 0.3 and e 2 = 0.7) (g2 was set to 0°). 
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Fig. 3. Transit timing variatons caused by a hypothetical P 2 = 10000 day-period = IMq mass, moderately eccentric e 2 = 0.3 third companion 
for CoRoT-9b, at different initial orbital elements, for four different initial inner eccentricities (e t = 0, 0. 1 1, 0.5, 0.9). The various initial elements 
of each panel are as follows. Panel (a): gi = 7°, g 2 = 90°, i m = 0°; (b): The same, but for i m = 90°; (c): g] = 7°, g 2 = 45°, i m = 30°; (d): gi - 337°, 
gi = 45°, ! m = 60°; (e): g\ = 307°, g 2 = 45°, i m = 30°; (/): same as previous, but for i m = 90°. (g): g\ = 277°, g 2 = 45°, i m = 0°; (h): same as 
previous, but for i m = 90°. (For better comparison the curves are corrected for the different average transit periods, and zero point shifts.) 
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where up to third order in e2 



Av2 ~ AZ2 ( 1 + 2e\ + 2e2 cos V2 + ^ e 2 cos 2v"2 + 3ej cos 3v2 



2tt—(1 +6v 2 ) 
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(54) 



i. e. we approximated the variation of the mean anomaly of the 
outer companion (fe) by the constant value of 



AZ 2 ~ 2tt 



Pi 



(55) 



(Here, and in the following text we neglect the pure geometri- 
cal LITE contribution.) Then the variation of the length of the 
consecutive transiting periods becomes 
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(56) 



Comparing (the amplitude of) this result with Eqs. (1,2) of 
Holman & Murray (2005), one can see that the power of the 
period ratio differs. (Furthermore, in the original paper n stands 
at the same place, i.e. in the numerator (as above) but later, in 
Holman 2010 this was declared as an error, and it was put into 
the denominator.) There is a principal difference in the back- 
ground. Holman & Murray (2005) state that they "estimate the 
variation in transit intervals between successive transits", but 
what they really calculate is the departure of a transiting period 
(i.e. the interval between two successive transits) from a constant 
mean value. This quantity was estimated in our Eq. (53). In other 
words, the Lagrangian perturbation equations gives the instanta- 
neous period of the perturbed system. So, its integration gives 
the time left between two consecutive minima. The first deriva- 
tive of the perturbation equation gives the instantaneous period 
variation, and so the variation of the length of two consecutive 
transiting period can be deduced by integrating the latter. As a 
consequence, the best possibility for the detection of the tran- 
sit timing variation, and so, for the presence of some perturber 
occurs when the absolute value of the second derivative of the 
O - C diagram, or practically Eq. (56) (the period variation dur- 
ing a revolution) is maximum. We will illustrate this statement 
in the next Section for specific systems. 



4. Case studies 

4.1. CoRoT-9b 

CoRoT-9b is a transiting giant planet which revolves around its 
host star approximately at the distance of Mercury (Deeg et al. 
2010). Consequently, tidal forces (including the possible ro- 
tational oblateness) can safely be neglected in this system. 
Furthermore, due to the relatively large absolute separation of 
the planet from its star, we can expect a large amplitude sig- 
nal from the perturbations of a hypothetical, distant (but not too 
distant) further companion. To illustrate this possibility, and, fur- 
thermore, to check our formulae, we both calculated and plotted 
the amplitudes with the measured parameters of this specific sys- 
tem, and carried out short-term numeric integrations for compar- 
ison. The physical, and some of the orbital elements of C0R0T- 
9b were taken from Deeg et al. (2010). These data are listed in 
Table 1. (Note, we added 180°to the u>\ published in that paper, 
as the spectroscopic oj\ refers to the orbit of the host star around 
the common centre of mass of the star-planet double system, and 
consequently, differ by 180°from the observational argument of 
periastron of the relative orbit of the transiting planet around its 
host star.) In the case of the numerical integrations, the inner 
planets were started from periastron, while the outer one from 
its apastron. 

Fixing the observationally aquired data, Am,s depend only 
on two parameters, namely gi and z' m . In the left panels of Fig. 4 
we plotted the i m versus Am ,s graphs for gi =69°andgi = 158°. 
We found that the amplitudes reach their extrema around these 
periastron arguments, i.e. for other gi values results occur within 
the areas limited by these lines. 6 In the middle and right panels 
the effect of the two further free parameters, i.e. outer eccen- 
tricity (e2), and dynamical (relative) argument of periastron (#2) 
were also considered. The middle panels show A 1,2 for e2 = 0.3, 
the right panels for e2 = 0.7, for both g2 = 0°(upper panels), and 
g2 — 90°(bottom panels). (Note, that As, plotted in the left panel 
is identical to A2 for e2 = 0, in which case A 13 = 0.) 

These figures clearly show again that the amplitudes, and 
consequently, the actual full amplitude of the dynamical O - C 
curve remains within one order of magnitude over a wide range 
of orbital parameters. This once more verifies the very crude es- 
timations given by Eqs. (51,52). As a consequence, for a given 
system we can estimate very easily the expected amplitude of the 
O - C variations induced by a further companion. For a planet- 
mass third body, the 



idyn 



1 m3 



2n P2 Whost-star 



(57) 



gives a likely estimation at least in magnitude. For example, for 
CoRoT-9b 



1 P 2 

m = ! — * HSgd^Q 1 

2n 

^host-star 

~ 1.39d 2 M J " 1 , 



(58) 



i.e. a Jupiter-mass additional planet could produce O?001 half- 
amplitude already from the distance of the Mars. (Of course, if 
e\ and co\ are known, a more precise estimation can be provided 
easily using the formulae of the present paper.) 

Nevertheless, while the net amplitudes usually vary in a nar- 
row range, the dominances of the two main terms (with period 



6 Strictly speaking this latter argument is true only if negative 
values are also allowed for the coefficients. 
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Fig. 4. Left panels: The mutual inclination (i ra ) dependence of Am, A$ amplitudes of M, S functions of the long-period dynamical component of 
O - C for CoRoT-9b, for g\ = 69°and g\ = 158°. For other gi values the corresponding curves run within the area limited by these lines. (The two 
left panels are identical.) Middle and right panels: The corresponding A l2 amplitudes of the trigonometric representation (Eq. 32) of the O-C 
for two different outer eccentricities e 2 = 0.3 (middle) and e 2 = 0.7 (right); and outer dynamical (relative) argument of pericentrum arguments 
g2 = 0°(up) and g 2 = 90°(down). For the sake of clarity, we did not plot the small A 3 coefficients. (Note, for e 2 = 0.0 the As amplitudes of the two 
identical left panels are equal to A?, while A] = A3 = 0.) 



P2 and can alternate and, although it is not shown, their 
relative phase can also vary so, the shape of the actual curves 
may show great diversity, as is illustrated e.g. in Figs. 3, 5. 
Furthermore, for specific values of the parameters, one or the 
other amplitudes might disappear. At CoRoT-9b particularly in- 
teresting is the g\ = 69°/ m ~ 45°configuration since in this case 
both Am and As disappear very close to each other. This means 
that for specific e^, gi values the dynamical O - C almost disap- 
pears. This possibility warns us of the fact, that from the absence 
of TTV in a given system one cannot automatically exclude the 
presence of a further planet (which should be observed accord- 
ing to its parameters). 

This can be seen clearly in Fig. 5, where we plotted the corre- 
sponding O-C curves for coplanar, perpendicular, and the above 
mentioned interesting i m = 40°and 46°configurations. In this lat- 
ter figure we plotted O-C curves obtained from both numerical 
integrations of the three-body motions, and analytical calcula- 
tions with our sixth order formula. (Note, according to Fig. 2, 
for the small inner eccentricity of CoRoT-9b (<?i = 0.11), the 
first order approximation would have given practically the same 
results.) Comparing the analytical and the numerical curves, 
the best similarity can be seen in the perpendicular case (last 
row). In the coplanar case some minor discrepancies can be ob- 
served both in the shape and amplitude, while the discrepancies 
are more expressed in the particular i m = 46°case, where for 
large outer eccentricity (right panel in the third row) our so- 
lution fails. (Nevertheless, the total amplitude of the analyti- 
cal curve is similar to the numerical one in this case, too.) 
Although a thorough analysis of the sources of the discrepancies 
is beyond the scope of this paper, we suppose that in these sit- 
uations the discrepancies come from the higher order perturba- 
tive terms. As was shown e.g. by Soderhjelm (1982); Ford et al. 



(2000) the higher order contributions are the most significant for 
ej ~ 0, i m = (0, 1) x 180°. Although the above authors con- 
sidered only the secular, or apse-node term perturbations, the 
same might be the case for the long-period ones. This might ex- 
plain the better correspondance in the perpendicular configura- 
tion than in the coplanar one. We suppose some similar reason 
in the i m = 46°case. As in this situation the first order contri- 
butions almost disappear, whereas the small higher order terms 
can also be more significant. Furthermore, numeric integrations 
in this latter case show a disappearence of the identity between 
the g\,gi + 180°initial conditions, which also suggests a signifi- 
cance of the higher order terms as in this case we can expect the 
appearance of trigonometric functions with g\ and 3g\ in their 
arguments. 

As expected, the highly eccentric-distant-companion sce- 
nario produces the largest amplitude TTV, at least when a com- 
plete revolution is considered. Nevertheless, on a shorter time- 
scale, the length of the observing window necessary for the de- 
tection depends highly on the phase of the curve. To illustrate 
this, in Fig. 6 we plotted the first, and second 8 years of the three 
primary transit O-C curves for both the coplanar, and the per- 
pendicular cases, shown in the first and last rows of Fig. 5. The 
transiting periods for each curve were calculated in the usual, ob- 
servational manner, i.e. the time interval between the first (some) 
transits were used. According to Fig. 6, in the present situation, 
in the first 8 years (i.e. which begins at the apastron of the outer 
planet) it would be unlikely to detect the perturbations in the 
largest (total) amplitude = 0.7 cases, as well as in the copla- 
nar <?2 = 0.3 case. The most certain detection would be possible 
in the two smallest amplitude circular — configurations. 
Nevertheless, if the observations starts at those phases plotted 
on the right panels (8-16 years), the pictures completely differ. 
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Fig. 5. Sample of transit timing variations caused by a hypothetical P% = 10000 day-period mj = 0.005 Mq(x 5Mj) mass third companion for 
CoRoT-9b at four different initial mutual inclinations (7 m = 0°, 40°, 46°, 90°, from up to down), and three different initial outer eccentricities 
(e2 = 0, 0.3, 0.7, from left to right). The dynamical (relative) arguments of periastrons are set to gi = 69°, g2 = 90°. The curves show the sum 
of the geometrical LITE, and the dynamical terms obtained both from numeric integrations, and analytic calculations up to sixth order in e\ . The 
pure LITE contributions are also plotted separately. Note that the vertical scale of the individual columns (i.e. different e2-s) are different. 



In this latter interval, the circular cases produce the smallest cur- 
vature O-C-s, and the discrepancy from the linear trend reaches 
$001 days (which can be considered as a limit for certain detec- 
tion) occuring towards the end of the interval. On the other hand, 
in the case of the highly eccentric configurations, the 'moment 
of truth' comes after some years. Nevertheless, we have to stress, 
that although we plotted the O - C curves with continuous lines, 
in reality they would contain only 3-4 points during the phase of 



the seemingly abrupt jump, is a further complicating factor with 
regards to 'certain' detection. 

In the sample runs above a moderately hierarchic sce- 
nario a 2 ~ 22a i was studied. In order to get some picture 
about the lower limit of the validity of our low-order, hi- 
erarchical approximation, we carried out further integra- 
tions for less hierarchic configurations. In Fig. 7 we show 
the results of some of these runs, which were carried out 
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Fig. 6. The first and the second 8 years of O - C-s plotted in the first and the last rows of Fig. 5. The periods of the individual curves were set equal 
to the respective initial transiting periods. 



with the same initial conditions what were used in the mid- 
dle panel of the first and last rows (coplanar and perpendic- 
ular cases, respectively) of Fig. 5, but for P 2 = 1000 days (i.e. 
a 2 ~ 4.8a!), P 2 = 2000 days (a 2 ~ 7.6^), and P 2 ~ 952.7 
days (a 2 ~ 4.6«i), in which latter case the two planets or- 
bit in 1 : 10 mean motion resonance. In the left panels of 
Fig. 7 we plot 20-year-long intervals, while in the right ones 
century-long time-scales are shown. As one can see, on this 
latter time-scale, apse-node effects already reach or exceed 
the magnitude of the long period ones. This naturally arises 
from the fact that the typical time-scales of these latter high- 
amplitude perturbations are proportional to P\lP\, i.e. in 
the present cases these are ~ lOOx faster than in the pre- 
viously investigated case. Note, that for the sake of a better 
comparison, the analytical curves in the right panels were 
calculated with the inclusion of these apse-node terms, al- 
though the latter will be presented only in a forthcoming pa- 
per. Turning back to the 20-year-long integrations, one can 
see that the limit of the validity of the present approxima- 
tion does strongly depend on the mutual inclination. While 
for the z' m = 90° configurations the long period analytical re- 
sults are in remarkably good agreement with the numerical 
curves even for a 2 < 5a 1 (left panels of third and forth rows), 
for the coplanar (i m = 0°) case our approximation is clearly 
insufficient for such small a 2 ja\ ratios, and even for the dou- 
bled outer period case (i.e. a 2 w 8a 1) the amplitude of the 
analytical curve is highly underestimated. 



We also investigated the case of the 1:10 mean-motion 
resonance. Our results for the perpendicular case are plot- 
ted in the last row. In this case the numerical integration 
shows very high amplitude apse-node scale variations that 
does not occur in the analytical curve. In order to get a better 
comparison between the analytical and numerical long-term 
variations in this case, we removed the apse-node effect from 
the numerical curve, by the use of a quadratic term, i.e. the 
(blue) O - C curve was calculated in the form of 

O - C = c + ciE + c 2 E 2 , (59) 

where E is the cycle number. As one can see, this quadratic 
(blue) curve shows similar agreement with the analytical 
curve, which was found in the similar 02/01 ratio non- 
resonant case. Consequently, we can state that our long term 
formulae are capable to produce the same accuracy even 
around mean-motion resonances. 

Nevertheless, from these few arbitrary trial runs we can- 
not give general statements about the limits of our approxi- 
mations. A detailed discussionof this point is postponed to a 
forthcoming paper, when we include the apse-node time scale 
terms. 

While CoRoT-9b served as an illustration for the Transit 
Timing Variations in the low inner eccentricity case, our next 
sample exoplanet HD 80606b represents the extremely eccentric 
case. 
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Fig. 7. Checking the validity of hierarchical approximation for closer systems. In the first two rows the initial conditions were set to be the 
same as at the uppermost middle panel of Fig. 5, with the exception of P 2 = 1000 (first row) and 2000 days (second row), i.e. a 2 la\ x 4.8 and 
« 7.6, respectively, while the third and last rows have initial conditions similar to the middle panel of Fig. 5, with the exception of P 2 = 1000 
(third row) and 952.7 days (last row). This latter illustrates the case of a 1 : 10 mean-motion resonance. The left panels represent a 20-year- 
long time-scale, while the right ones show the TTV behaviour during a century. In the left panel of the last row the blue (quadratic) curve 
shows the O - C curve calculated by including a quadratic term. See text for details. 
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Table 1. The initial parameters of the transiting planetary subsytems. 
(The masses are given in solar mass, periods in days, and angular ele- 
ments in degrees.) The parameters are taken from Deeg et al. (2010) for 
CoRoT-9b, and from Pont et al. (2009) for HD 80606b. (Note that the 
argument of periastron (a>i) for the relative orbit of the planet around 
its host star differ by 180° from the value deduced from radial velocity 
data.) 



System 


mi 


m 2 


Pi 


e\ 


a>i 


h 


CoRoT - 9b 
HD 80606b 


0.99 
0.97 


0.0008 
0.0038 


95.2738 
111.4357 


0.11 
0.93 


217 
121 


89.99 
89.32 



4.2. HD 80606b 

The high-mass gas giant exoplanet HD 80606b features an al- 
most 4-month long period and an extremely eccentric orbit 
around its solar-type host-star. It was discovered spectroscopi- 
cally by Naef et al. (2001). Recently both secondary occultation 
(with the Spitzer space telescope, Laughlin et al. 2009), and pri- 
mary transit (Moutou et al. 2009; Fossey et al. 2009) have been 
detected. A thorough analysis of the collected data around the 
February 2009 primary transit led to the conclusion that there 
is a significant spin-orbit misalignment in the system, i.e. the 
orbital plane of HD 80606b fails to coincide with the equato- 
rial plane of its host star (Pont et al. 2009). These facts suggest 
that this planet might be seen at an instant close to the maxi- 
mum eccentricity phase of a Kozai cycle induced by a distant, 
inclined third companion (cf. Wu & Murray 2003; Fabry cky & 
Tremaine 2007). Note that although HD 80606 itself forms a bi- 
nary with HD 80607, due to the large separation, one may expect 
a further, not so distant companion for an effective Kozai mech- 
anism. Such a very distant companion may nevertheless play an 
indirect role in the initial triggering of the Kozai mechanism, by 
the effect described in Takeda et al. (2009). The most important 
parameters of the system are summarized in Table 1 . 

Due to its very high eccentricity, the periastron distance of 
this planet is qi = ai(l — ei)« 0.03 AU at which distance the 
tidal forces and especially tidal dissipation should be effective. 
(As a comparison we note, that this separation corresponds to the 
semi-major axis of an approx. 2 day period orbit.) Furthermore, 
we can expect also a significant relativistic contribution to the 
apsidal motion. Although such circumstances do not invalidate 
the following conclusions (as the time-scale of these effects are 
significantly longer than the ones investigated) we nevertheless 
provide a quantitive estimation of their effects. Furthermore, the 
effect of dissipation will also be considered numerically, at the 
end of this section. 

As is well known, the apsidal advance speed, averaged for 
one orbital revolution, can be written in the following form: 



gi = A + Bcos2gu 



(60) 



where the non-zero contributions of the third-body, tidal and rel- 
ativistic terms are as follows: 
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In these equations mi, Ri, pi, k^\ w z <, v e , refer to the mass, ra- 
dius, average density, first apsidal motion constant, uni-axial ro- 
tational angular velocity, and equatorial rotational velocity of the 
host star respectively, while subscript 2 denotes the same quanti- 
ties for the inner planet. Note also that the rotational term is valid 
only for non-aligned rotation, and consequently, it provides only 
a crude estimation for HD 80606b. Nevertheless, for the present 
purpose it seems satisfactory. First consider the third-body term. 
Substituting cos2g = -1, 1 2 = 3/5, 



(68) 
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from which the last term usually omittable in hierarchical sys- 
tems, since C\ « C 2 , which is more expressively true for high 
ei . For example, in the present situation 



So, for HD 80606b we obtain that: 
Agy* * 0?43 x (l - eiy V2 century" 1 . 



(70) 



(71) 



(This result is in excellent correspondance with Fig. 11.) 

There are several uncertainties in the calculation of the tidal 
contribution. While the k 2 constant is relatively well-known 
for ordinary stars, it has a great ambiguity for exoplanets. 
Furthermore, we do know nothing about the rotational velocity 
of HD 80606b. So, according to the tables of Claret & Gimenez 
(1992) we set = 0.02 for the host star, and assumed k^ = 0.2 
for its planet, which is the same order of magnitude as for Jupiter 
and of WASP-12b (Campo et al. 2010). The stellar rotation was 
set to V rot = 1.8 kins" 1 (i.e. P rot = 27<?5; oj z[ = 0.228 day" 1 ) 
(Fischer & Valenti 2005), while for the planet we supposed (ar- 
bitrarily) a one-day rotation period, i.e. ojj = 6.283 day -1 . By 
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the use of these values, the classical tidal contribution to the ap- 
sidal motion becomes 

Agddd ~ 0°007 century" 1 , (72) 

which can be neglected. 

Finally, the relativistic contribution is estimated to be: 

Ag re i ~ 0°06 century" 1 , (73) 

i.e. it is smaller by one magnitude than the third-body term, and 
consequently, does also not play any important role. 

Returning to the P 2 time-scale variations of the TTV-s, we 
carried out our calculations and integration runs with the suppo- 
sition that a second, similar mass giant planet is the source of 
this comet-like orbit, which is seen in the instant of the max- 
imum eccentricity phase of the Kozai-cycle. Consequently, we 
set gi = 90°, and i m = 39°23. With these values from the sixth 
order formula we get Am = -0.38, As = 1.69 for primary 
transits, and Am = -1.06, As = 1.72 for secondary occulta- 
tions. Consequently, for primary transits the O - C is evidently 
dominated by the <S-term. (The negative Am indicates a simple 
180°phase-shift.) In Fig. 8 the A 1,2,3 amplitudes are plotted as 
a function of the outer eccentricity (ei). Due to the more than 
four and a half-times larger primary transit S amplitude, the g2 
dependence of the amplitudes here are weak, and therefore we 
show only Ai^-s for g2 = 0°. In Fig. 9 we present both the 
numerically generated short-term O - C curves, as well as the 
analytically calculated cases up to the sixth order in eccentric- 
ity (see Appendix) for three different eccentricities (e 2 = 0, 0.3, 
0.7) of the outer perturber's orbit. We plotted the O-C-s for both 
primary transits and secondary occultations. Also shown are the 
corresponding primary minus secondary curves. 

As one can see, the sixth order formulae gave satisfactory 
results even for such high eccentricities, although the analytical 
amplitudes are somewhat overestimated. Nevertheless, a more 
detailed analysis shows that the accuracy of our formulae for 
such high inner eccentricities strongly depends on the other or- 
bital parameters. This is illustrated even in the present situation, 
where, for the secondary-occultation curves, (which correspond 
to a>i = 301°), the discrepancies are clearly larger. From this 
point of view, the e 2 = case (first row) is the more interesting, 
as in this situation, due to the non-zero, very similar As = Ai am- 
plitudes, we would expect almost identical O - C curves (since 
the dashed analytical ones are very similar), but this, in fact, is 
not the case. 

Considering the case of fastest possible detection of the am- 
plitudes, in Fig. 10 we plotted also the first and second 8 years of 
the three primary transit O-C curves, shown in Fig. 9. The tran- 
siting periods for each curves were calculated on the usual, ob- 
servational manner, i.e. the time interval between the first (some) 
transits were used. According to Fig. 10, in the first 8 year, i.e. 
after the apastron of the outer body, the fastest detection would 
be possible in the smallest (total) amplitude circular third body- 
case, while the curvature of the highly eccentric curve is so 
small, that it needs almost 7 years to exceed the 0?001 differ- 
ence which could promise certain detection. (Of course, this is a 
non-realistic ideal case, when all the transits are measured, with- 
out any observational error.) Around periastron (right panel), the 
situation is completely different, similarly to the case of C0R0T- 
9b. 

We also carried out numerical integrations to investigate the 
possible orbital evolution of HD 80606b in the presence of such 
a perturber. We have shown that both relativistic and tidal ef- 
fects are omittable in the present configuration of the system. 



However, this assumption is only correct in those cases where 
tidally forced dissipation is not considered. Consequently, we 
integrated the dynamical evolution of the system including both 
tidal effects (both dissipative and conservative tidal terms), and 
without them (i.e. in the frame of pure three, point-mass grav- 
itational interactions). The equations of the motion (including 
tidal and dissipative terms, and stellar rotation) are given in 
Borkovits et al. (2004), where the description of our integrator 
can also be found. In the dissipative case our dissipation constant 
was selected in such a way that it produced At\ ~ 2^16 x 10~ 7 
tidal lag-time for the host-star, and Afi « 5*?21 x 10~ 4 for the 
planet, which are equivalent to Q\ « 4.1 X 10 7 , Q 2 « 1.7 X 10 4 
dissipation parameters, respectively. Fig. 1 1 shows the variation 
of most of the orbital elements (both dynamical and observa- 
tional) for 100 and 1 million years in the dynamically most- 
excited e2 = 0.7 case. Thick lines represent the dissipative case, 
while thin curves show the point-mass result. The century-long 
left panel demonstrates, that apart from the shrinking semi-major 
axis, there is no detectable variation in the orbital elements dur- 
ing such a short time-interval. And, of course, if this is true for 
the extremum of the Kozai cycle, it is more expressly valid for 
other situations, as this phase produces the fastest orbital ele- 
ment variations. Furthermore, on such a short time-scale the or- 
bital variations in a point-mass or non-Keplerian, tidal frame- 
work are indistinguishable. (Again, we do not take into account 
the semi-major axis.) This provides further verification of the 
effects previously-discussed, namely that we neglected the tidal 
and relativistic effects completely, and considered all the orbital 
elements as constant. 

Now, we consider orbital shrinking due to dissipation. As 
one can see, in the present situation the decrease in a\ during the 
first 100 years is ha\ ~ 6x 10~ 4 R©. Converting this into a period 
variation suggests APi « 10~ 3 d. This gives for one transiting 
period a Pi « 3x 10~ 6 day cycle -1 rate. From the point of view of 
an eclipsing binary observer, this is an incredibly large value. For 
comparison, a typical secular period variation rate measured for 
many of the eclipsing binaries is about 10~ 9 - 10~ n day cycle" 1 . 
Such a high rate would produce 0?001 departure in transit time 
during ~ 26 cycles, i.e. during approximately 8 years. Note, this 
period variation ratio is close to that ratio produced by typical 
long-term perturbations within a few years. This illustrates that 
if the data-length is significantly shorter than the period of the 
long-term periodic perturbations, then the effect arose from such 
perturbations, and other effects, coming from i.e. orbital shrink- 
ing can overlap each-other, and as a result they can be misinter- 
preted easily. (The question of such kinds of misinterpretations 
or false identifications were considered in general in Sect. 2 of 
Borkovits et al. 2005). 

Finally, we consider the right panel of Fig. 11, which il- 
lustrates the frequently mentioned Kozai-mechanism (with and 
without tidal friction) in operation. Note, that in the present sit- 
uation due to the high outer eccentricity, and as well as the rel- 
atively weak hierarchicity of the system (i.e. a\/ci2 ~ 0.05) the 
higher order terms of the perturbation function are also signifi- 
cant, which results in very different consecutive cycles even in 
the point-mass case, too. This manifests not only as different pe- 
riods and maximum eccentricities, but e.g. in the fact that the 
argument of (dynamical) periastron (gi) shows both circulation, 
and libration, alternately (cf. Ford et al. 2000). Nevertheless, the 
detailed investigation of the apse-node timescale behaviour of 
the TTV, as well as other orbital elements and observables will 
be the subject of a succeeding paper. 
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Fig. 8. The A 1.2,3 amplitudes for HD 80606b for primary transits (left), and secondary occultations (right), supposing that the planet is seen at the 
instant of the maximum phase of a Kozai-cycle. 



5. Conclusions 

We have studied the long-term P2 time-scale transit timing vari- 
ations in transiting exoplanetary systems which feature a further, 
more distant (ai » a\) either planetary, or stellar companion. 
We gave the analytical form of the O - C diagram which de- 
scribes such TTV-s. Our result is an extention of our previous 
work, namely Borkovits et al. (2003) for arbitrary orbital ele- 
ments of both the inner transiting planet and the outer compan- 
ion. We showed that the dependence of the O - C on the or- 
bital and physical parameters can be separated into three parts. 
Two of these are independent of the real physical parameters (i.e. 
masses, separations, periods) of a concrete system, and depend 
only on dimensionless orbital elements, and so, can be analysed 
in general. For the two other kinds of parameters, which are am- 
plitudes and phases of trigonometric functions, we separated the 
orbital elements of the inner from the outer planet. The practical 
importance of such a separation is that in the case of any ac- 
tual transiting exoplanets, if eccentricity (<?i) and the observable 
argument of periastron {lo\) are known e.g. from spectroscopy, 
then the main characteristics of any, caused by a possible third- 
body, transit timing variations can be mapped simply by the vari- 
ation of two free parameters (dynamical, relative argument of 
periastron, g\, and mutual inclination, i m ), which then can be 
refined by the use of the other, derived parameters, including 
two additional parameters fot the possible third body (eccen- 
tricity, ez, and dynamical, relative argument of periastron, g2). 
Moreover, as the physical attributes of a given system occur only 
as scaling parameters, the real amplitude of the O — C can also be 
estimated for a given system, simply as a function of the m^jPz 
ratio. 

At this point it would be no without benefit to compare 
our results with the conclusions of Nesvorny & Morbidelli 
(2008) and Nesvorny (2009). These authors investigated the 
same problem, i.e. the fast detectability of outer perturbing 
planets, and determination of their orbital and physical pa- 
rameters from their perturbations on the transit timing of 
the inner planet, by the help of the analytical description of 
the perturbed transiting O - C curve. For the mathemati- 
cal description they used the explicite perturbation theory of 
Hori (1966) and Deprit (1969) based on canonical transfor- 
mations and on the use of Lie-series. This theory does not 
require the hierarchical assumption, i.e. the a - ai/a 2 pa- 
rameter, although less than unity, is not required to be small. 



This suggests a somewhat greater generality of the results, 
i.e. it is well applicable systems similar to our solar system. 
Nevertheless, perhaps a small disadvantage of this method 
with respect to the hierarchical approximation is, that for 
large mutual inclinations the number of terms in the per- 
turbation function needed for a given accuracy grows very 
fast, while our formulae have the same accuracy even for the 
largest mutual inclinations. There is also a similar discom- 
fort in the case of high eccentricities, in which case the gen- 
eral formulae of Nesvorny (2009) are more sensitive for these 
parameters than in the hierarchical case. For example, we 
could reproduce satisfactory accuracy even for e\ = 0.9, in 
which case the classical formulae are divergent. (Note, that 
although in this paper we concentrate only on the long pe- 
riod perturbations, the same is valid for the apse-node time- 
scale variations, as will be illustrated in the next paper.) As 
a conclusion, the greater generality of the method used by 
Nesvorny & Morbidelli (2008) and Nesvorny (2009) beyond 
the evident non-hierarchical configurations is well (or even 
better) applicable also in the nearly coplanar case, especially 
when the inner orbit is nearly circular (in which strict case 
the first order hierarchical approximation becomes insuffi- 
cient, see e.g Ford et al. 2000), but in other cases, as far as 
the hierarchical assumption is satisfied, this latter could give 
a faster and simpler method. Furthermore, since in the hier- 
archical approximation, the principal small parameter in the 
perturbation equations is the ratio of the separations instead 
of the mass-ratio, these formulae are valid for stellar mass 
objects as well, and can also be applicable for planets orbit- 
ing an S-type orbit in binary stars, or even for hierarchical 
triple stellar systems. 

We analysed the above-mentioned dimensionless amplitudes 
for different arbitrary initial parameters, as well as for two con- 
crete systems CoRoT-9b and HD 80606b. We found in general, 
that while the shape of the O - C strongly varies with the an- 
gular orbital elements, the net amplitude (departing from some 
specific configurations) depends only weakly on these elements, 
but strongly on the eccentricities. Nevertheless, we found some 
situations around i m = 45°for the specific case of CoRoT-9b, 
where the O - C almost disappeared. 

We used CoRoT-9b and HD 80606b for case studies. Both 
giant planets revolve on several month-period orbits. The former 
has an almost circular orbit, while the latter has a comet-like, ex- 
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Fig. 9. Transit timing variatons caused by a hypothetical P 2 = 10000 day-period 1113 = 0.005Mq(»; 5Mj) mass third companion for HD 80606b at 
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tremely eccentric orbit. These large period systems are ideal for 
searching for further perturbing components, as the amplitude 
of the O - C is multiplied by P\jP2 and consequently, as the 
magnitude of the perturbations determined by Pi/Pz, the same 
amplitude perturbations cause a better detectable effect, if the 
characteristic size of the system (i.e. Pi) is larger. 

We considered also the question of detection, as well as the 
correct identification of such perturbations. We emphasize again, 
that the O-C curve is a very effective tool for detection of 
any period variations, due to its cumulative nature. Nevertheless, 
some care is necessary. First, it has to be kept in mind, that the 
detectability of a period variation depends on the curvature of 
the O-C curve. (If the plotted O - C is simply a straight line, 
with any slope, it means that the period is constant, which is 
known with an error equal to that slope.) This implies that the 
interval which is necessary to detect the period variation com- 
ing from a periodic phenomenon depends more strong on which 
phase is observed than on the amplitude of the total variation. 
(We illustrated this possibility for both system sampled.) We il- 
lustrated also, that in the case of a very eccentric third compan- 
ion the fastest rate period variation lasts a very short interval, 
which consists of only a few transit events. This emphasizes the 
importance of observing all possible transits (and, of course oc- 



cultation) events, with great accuracy. A further question is the 
possible misinterpretation of the O-C diagram. For example, 
as was shown, in the HD 80606b system we can expect a sec- 
ular period change due to tidal dissipation, which has the same 
order of magnitude than might have been measured due to the 
periodic perturbations of some hypothetical third body as our 
sample. These two types of perturbations could be separated in 
two different ways. One way is simply a question of time. In the 
case of a sufficiently long observing window, a periodic pertu- 
bation would separate from a secular one, which would produce 
a parabola-like O-C continuously. But, there is also a faster 
possibility. In the case that not only primary transits, but also 
secondary occultations are observed with similar frequency and 
accuracy, then, on subtracting the two (primary and secondary) 
O-C curves from each-other, the secular change, i.e. the effect of 
dissipation would disappear, since it is similar for primary tran- 
sits and secondary occultations. This fact also makes desirable 
the collection of as many as possible occultation observations 
too. 
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Appendix A: Relation between the observable, and 
the dynamical orbital elements. 

The start of the paper remarked that some of the relations be- 
tween the elements given in this section are valid strictly in the 
presented form only for the situation shown in Fig. 1 . According 
to the actual orientations of the orbital planes, the spherical trian- 
gle (with sides u m \, u m 2, Qi - Q2) could be oriented in different 
ways, and so, some precise discussion are necessary, which is 
omitted in the present appendix. 

The relation between the pericentrum arguments are as fol- 
lows: 

«1 = gl + "ml, (A.l) 
«2 = g 2 + Km2+180°. (A.2) 

Consequently, the true longitudes measured from the sky (m) and 
from the intersection of the two orbital planes (w) are: 



"1 


— Vi + COi, 


(A.3) 


"2 


- V 2 + Cl>2, 


(A.4) 


W] 


= Vl +gu 


(A.5) 




= U\ — U m i, 


(A.6) 


\V 2 


= v 2 +g 2 + 180°, 


(A.7) 




- U2~ W m 2- 


(A.8) 



There are two relations between the mutual inclination (i m ) and 
the two dynamical inclinations (Ji,2)- One of these is trivial, 
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Fig. 11. Dynamical evolution of the orbital elements of HD 80606b in the presence of a hypothetical P 2 = 10 000 day-period m 3 = 0.005Mq(«; 
5Mj) mass third companion. The initial orbital elements correspond to the last row of Fig. 9. Red curves: tidal effects and dissipation are considered; 
blue curves: three point-mass model. (Relativistic contributions are omitted.) Note, the semi-major axis (a t ) is given in Rq. 



while the second comes from the relation of the inclinations to 
the orbital angular momenta. The trivial case is 



*m = J\ + J2, 

while the non-trivial case comes from the fact, that 
CCi 



COS J] = 

sin ji = 
i.e. 

COSJl = 

sin ./'i = 



ccy 

ICxCil 
CCi ' 

Ci c 2j 
-C+-C 1 ' 



c 



■ sin i n 



(A.9) 

(A. 10) 
(A.ll) 

(A. 12) 
(A.13) 



where represents the orbital angular momentum vector of 
the two orbits, C the net orbital angular momentum vector, and 



C, 



c 2 



m\m2 



^Gm n a y (l -efj, 
- s jGmn3a2 (l - 4), 



mm 

C = C\ cos j\ + C2 cos 7*2 



(A. 14) 

(A.15) 
(A. 16) 



(The rotational angular momenta were neglected.) At this point 
we remark that for hierarchical systems it is a very reasonable to 
assume 



Ci cos j\ - constant, 



(A. 17) 



(since the first order, doubly averaged Hamiltonian of such sys- 
tems does not contain its conjugated variable, h\), which prop- 
erty, together with the constancy of the semi-major axis connects 
the eccentricity (ei) with the dynamical inclination (/1). 

Further relations can be written by the use of the several 
identities of spherical triangles. For example: 

cos i m = cos i\ cos z'2 + sin z'i sin 12 cos(Q2 - Oi), (A. 18) 
sin i m sin M m2 = sin i\ sin(Q2 - £^i), (A. 19) 

sin/ m sinM m i = sin z'2 sin(D 2 - ^1), (A. 20) 

sin z m cos M m 2 = - cos z'i sin/ 2 + sin z'i cos z 2 cos(Q2 - QiXA.21) 
sin i m cos u m \ = cos z 2 sin z'i - sin z 2 cos i\ cos(Q 2 - Qi), (A. 22) 

and similar ones can be written for the two smaller spherical 
triangles. Finally, in hierarchical systems usually C 2 » C\, so 
we can take the following as a first approximation: 



71 

h 



(A.23) 
(A.24) 



22 



T. Borkovits et al.: Transit timing variations in eccentric hierarchical triple exoplanetary systems 



'0 - '2, 
hi = u m2 . 



(A.25) 
(A.26) 



We note that in Fig. 1, for practical reasons, we chose the 
intersection of the invariable plane and the sky for the arbitrary 
starting point of the observable nodes (Oi, Q2). Traditionally, in 
astrometry, these quantities are measured from north to east on 
the sky, nevertheless such an approximation is valid since the Q- 
s only appear in the equations in the form of their differences, 
and derivatives. (Nodes cannot be determined from photometry 
and spectroscopy, as both the light curves and the radial veloc- 
ity curves are invariant for any orientation of the orbital planes 
projected onto the sky.) 
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Appendix B: The derivatives used for calculating the long-period dynamical O - C up to sixth order in e\ 



The derivatives of the indirect terms are as follows: 
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|/ 2 -iJ + (l-/ 2 )cos(2v 2 + 2g 2 ) 

+6e^cos(6wi - 2gi) [(l - / 2 ) + (l + 7 2 )cos(2v 2 + 2g 2 )] 
-6e\ sin(6wi - 2gr)27sin(2v 2 + 2g 2 )} 

+ sim m cot n f + 3 2j COSMmi + 6e 2 cos(2gi + M m i) e{ cos 6a>il [1 - cos(2v 2 + 2g 2 j] 

(l-e 2 )' U 5 \ 2 / J 



12 



— 1 + -e\\ sin w m i + 6e 2 sin(2gi + w m i) 



e\ cos 6wi sin(2v2 + 2g 2 ) 



}}■ 



5 \ 2 

The direct terms are coming from (a 3 ^ 2 e m cos(nv) and Cl cos /1) are as follows: 

(^F) = A L (l+ e2 cosv 2 )(l- e 2 ) 1/2 |-^l + i e 2 J |/ 2 -iJ + (l-/ 2 )cos(2v 2 +2g2) 

-e 2 cos 2gi [(l - / 2 ) + (l + 7 2 ) cos(2v 2 + 2g 2 )] 
-e\ sin 2gi27 sin(2v 2 + 2g 2 )} + 



Me 2 / 3 (ei)cos2v '| 
V dv 2 J dir 



/ dej/6(ei)cos3v 

dV2 



du- 



de j/g(ei) cos 4 v 

dV2 



dir 



A L (1 + e 2 cos v 2 ) (1 - e\) 111 \^e\ |l + i e 2 + i^?J (/ 2 - I J + (l - 7 2 ) cos(2v 2 + 2g 2 ) 

-\ e2 \ l - % e \ - l e i) cos 2 Si [(l - I 2 ) + (l + cos(2v 2 + 2g 2 )] 
-^e 2 |l - ^e 2 - ^|Jsin2^ 1 2/sin(2v 2 + 2g 2 )J + 

A L (l+e 2 cos v 2 ) (1 - e\) 112 l- 3 -e\ (l + i e 2 J (/ 2 - I J + (l - I 2 ) cos(2v 2 + 2g 2 ) 
+ ^e\[\+e\- ^e A \ cos 2 gl [(l - 7 2 ) + (l + 7 2 ) cos(2v 2 + 2g 2 )\ 
+ ^e 2 {l +e\- ^etjsin2 <?1 2/sin(2v 2 + 2g 2 )j + -, 

A L (1 + e 2 cos v 2 ) (1 - e\) 112 |Ie« |/ 2 - I J + (l - 7 2 ) cos(2v 2 + 2g 2 ) 
-|ef |l + ^ 2 Jcos2g, [(l - 7 2 ) + (l + 7 2 )cos(2v 2 + 2* 2 )] 



"I** I 1 + ^'il s ^i2/sm(2r 2 + 2g 2 )\ + .... 



/ dej cos5v 



dv 2 



dir 



= A L (1 +e 2 cos v 2 )(l -e 2 ) 1/2 |ie«cos2gi [(l -7 2 ) + (l + 7 2 ) cos(2v 2 + 2# 2 )] 



+ -e\ sin2gi27sin(2v 2 +2g 2 )\ + 
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where +... refer to those terms which come from the normal force component, and will be cancelled by equal but opposite in sign 
direct Cl- terms. Furthermore, from 



3/2 \ 



dv2 (1 + e\ cos v) 2 



'-^A L (1 + e 2 cos v 2 ) (1 - ef) 1/2 \^e\ 1 1 + — e\ + —e\ 



2n 



/dir 



9 2 , 7 . 



16 



16 



/ 2 - - +(l -/ 2 )cos(2v 2 +2g 2 ) 



+ § e ' f 1 + t/' + TE eA ) cos 281 K 1 " /2 ) + ( : + /2 ) C0S(2V2 + 2g2) J 



21 



20 



57 



14 



112 



fjsir 



+ ™ e T 1 + T7 e i + TP7 g i sin 2^i2/sin(2v 2 + 2g 2 ) 



and from £l-term: 



(1-4) 



3/2 x 



dO 

COSii 

dv 2 (1 + e\ cos v) z 



A L (1 + e 2 cosv 2 ) 



/dir 



sin / m cot ;'i 

d-e 2 ) 1/2 



1 + -ef cos M mi + ef cos(M m i + 2gi) 



-- 1 + -ef sin M m , + ef sin(t< m i + 2g\) 



By the use of equations above, and integrating for v 2 , finally we get: 



sin(2v 2 + 2g 2 ) 



}- 



(B.13) 

I [I - cos(2v 2 + 2g 2 )] 
(B.14) 



(0-Q V2 = ^ L <|(l-e 2 ) 



2tt 



,1/2 



4 6 

1 



(/ 2 -^M+^(l-/ 2 )^(2v 2 + 2^ 2 ) 



51 2 

^<?i/20i)cos2,gi + 2/T 2 (e 1 ,wi,,gi) + -e^ieu^ug^ 



51 , 1 , 

^ei/2(ei)sin2gi + 2K 3 (e u co u gi) + -e\K 5 {e x ,u) U g{) 



sin; m cotii 

Mi* 

2 
5 



- 1 + -ej cos M ml + e l cos(2gi + w m i) 



(l-/ 2 )M+i(l+/ 2 )^(2v 2 + 2^ 2 ) 
2/C(2v 2 + 2g 2 )J 

[1+2^(^,^)1/ 



1 



M - -S(2v 2 + 2g 2 ) 



1 

+ - 

2 



sinM m i +e!sin(2gi +u m \) 



[\+2K l (e u u l )]C{2v 2 + 2g 2 ) 



where 



3 1 5 3-7 

K\{ei,a>\) = +e\ sin^! + -ef/3 cos2<di ± — g]/ 6 sin3wi - y^ e i/8 cos4wi + — ef sin5wi + — ef cos6wi, 

3 15 

K 2 (ei,ui,gi) = +eisin(wi -2gi)+ -e\f A cos(2wi -2g\)± -ef/ 7 sin(3wi - 2gi) - — ef/ 5 cos(4wi - 2gi) 

4 z lo 

3 7 
+— ef sin(5wi - 2gi) + — ef cos(6w! - 2gi), 
16 64 

3 15 

£3(21. wi.gi) = +e\ cos(w! - 2gi) - -e\f A sin(2w! - 2gi) ± -ef/ 7 cos(3wi - 2gi) + — ef/ 5 sin(4wi - 2g0 

4 z lo 

3 - 7 
+— cos(5w! - 2gi) - — ef sin(6wi - 2g\), 

3 

# 4 (ei,wi,gi) = -e 2 / 5 cos(2wi + 2g{) + ef sin(3^i + 2gi) + -efcos(4wi + 2g{), 

3 

/: 5 (ei, = -e 2 / 5 sin(2wi + 2g x ) ± ef cos(3^i + 2gi) + -ef sin(4wi + 2gi), 

and 

. , 25 2 15 4 95 6 
fi = 1 + — ef + — e, + — e, , 
71 8 1 8 1 64 1 



51 



48 1 



/3 = 1 + r 2 + ^>' 

1 2 1 4 

/ 4 = 1 + -e x + -e 1? 



(B.15) 

(B.16) 

(B.17) 

(B.18) 
(B.19) 
(B.20) 

(B.21) 
(B.22) 
(B.23) 
(B.24) 
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fs = 1 + \e\, (B.25) 

h = 1 + (B.26) 

fi = 1 + (B.27) 

/ 8 = 1 + (B.28) 
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